/******************************************************************************\
 *                                                                            *
 *                          Craftmesh model                                   *
 *                   Flavio Bertini fl.bertini@gmail.com                      *
 *                                                                            *
 *                                                                            *
\******************************************************************************/


#include "Craftmesh_model.h"
#include "Craftmesh_gengine.h"

using namespace Leap;
using namespace std;


// Create a mesh
vector<Mod_SubdivMesh> _mesh;


// Create an empty model
Mod_SubdivMesh::Mod_SubdivMesh()
{
  // init
  this->filename.append("");
  this->filepath.append("");
  this->nVert = 0;
  this->nFace = 0;
  this->vert  = NULL;
  this->face  = NULL;
  this->fNorm = NULL;
  mat_glMatrixfIdentity(this->scale);
  mat_glMatrixfIdentity(this->invsca);
  mat_glMatrixfIdentity(this->rotat);
  mat_glMatrixfIdentity(this->invrot);
  this->color = GL_COLOR_V[0];
  this->handId = NULLREF;
  this->captured = false;
  return;
}

// Create a model
Mod_SubdivMesh::Mod_SubdivMesh(string fName, string fPath)
{
  // init
  this->filename.append(fName);
  this->filepath.append(fPath);
  this->nVert = 0;
  this->nFace = 0;
  this->vert  = NULL;
  this->face  = NULL;
  this->fNorm = NULL;
  mat_glMatrixfIdentity(this->scale);
  mat_glMatrixfIdentity(this->invsca);
  mat_glMatrixfIdentity(this->rotat);
  mat_glMatrixfIdentity(this->invrot);
  this->color = GL_COLOR_V[0];
  this->handId = NULLREF;
  this->captured = false;
  if (!mod_load(fName, fPath))
  {
    cout << "Error: Unable to load " << filename << " -> exit!!" << endl;
    exit(1);
  }
  return;
}

// Destructor
Mod_SubdivMesh::~Mod_SubdivMesh()
{
  /* Se si libera la memoria qui, quando si aggiungono oggetti mesh a un vector ci sono problemi.
   * Infatti, se il vector fa un resize in seguito ad una push_back, gli oggetti vengono copiati
   * nel nuovo vector e i vecchi distrutti. Ma le nuove copie hanno gli stessi puntatori dei vecchi
   * oggetti: non si può liberare la memoria con il distruttore!!! */
  // free(this->vert);
  // free(this->face);
  // free(this->fNorm);

  return;
}

// load the obj model
bool Mod_SubdivMesh::mod_load()
{
  return mod_load(filename, filepath);
}

bool Mod_SubdivMesh::mod_load(string fName, string fPath)
{
  // obj wavefront to load
  this->filename.append(fName);
  this->filepath.append(fPath);
  string file = "";
  file += fPath;
  file += fName;
  file += ".obj";

  // file descriptor opening
  ifstream fd;
  fd.open(file.c_str());
  if (!fd.good())
  {
    cout << "Error: Unable to open " << filename << " -> exit!!" << endl;
    exit(1);
  }
  fd.seekg(0, ios::end);
  if (!fd.tellg())
  {
    cout << "Error: The file " << filename << " is empty -> exit!!" << endl;
    exit(1);
  }
  fd.seekg(0, ios::beg);

  // model statistics: num vertices and num faces
  string line = "";
  while (!fd.eof())
  {
    getline(fd, line);
    if (line.c_str()[0] == 'v') nVert++;    // load num vertices
    if (line.c_str()[0] == 'f') nFace++;    // load num faces
  }
  fd.close();

  // model elements
  fd.open(file.c_str());
  vert = (Vertex*) malloc((nVert+1)*sizeof(Vertex));
  if (!vert)
  {
    cout << "Error: Malloc fails -> exit!!" << endl;
    exit(2);
  }
  face = (Face*) malloc((nFace+1)*sizeof(Face));
  if (!face)
  {
    cout << "Error: Malloc fails -> exit!!" << endl;
    exit(2);
  }
  fNorm = (Vertex*) malloc((nFace+1)*sizeof(Vertex));
  if (!fNorm)
  {
    cout << "Error: Malloc fails -> exit!!" << endl;
    exit(2);
  }
  while (!fd.eof())
  {
    getline(fd, line);

    // load vertices
    if (line.c_str()[0] == 'v' && line.c_str()[1] == ' ')
    {
      for (int i=0; i<nVert; i++)
      {
        if (i) getline(fd, line);
        line[0] = ' ';
        sscanf(line.c_str(), "%f %f %f ", &(vert[i+1].x), &(vert[i+1].y), &(vert[i+1].z));
      }
      mod_setCenterAndRadius();
    }
    // load faces
    if (line.c_str()[0] == 'f')
    {
      for (int i=0; i<nFace; i++)
      {
        if (i) getline(fd, line);
        vector<int> buffer;
        int vIndex  = 0;
        int valency = 0;
        line[0] = ' ';
        stringstream indexList(line);
        while (!(indexList >> vIndex).fail())
        {
          buffer.push_back(vIndex);
          valency++;
        }
        face[i+1].nVert = valency;
        for (int j=0; j<valency; j++)
          face[i+1].vert[j] = buffer[j];
        // compute face normal
        float triangle[9];
        triangle[0] = vert[face[i+1].vert[0]/*-1*/].x;
        triangle[1] = vert[face[i+1].vert[0]/*-1*/].y;
        triangle[2] = vert[face[i+1].vert[0]/*-1*/].z;
        triangle[3] = vert[face[i+1].vert[2]/*-1*/].x;
        triangle[4] = vert[face[i+1].vert[2]/*-1*/].y;
        triangle[5] = vert[face[i+1].vert[2]/*-1*/].z;
        triangle[6] = vert[face[i+1].vert[1]/*-1*/].x;
        triangle[7] = vert[face[i+1].vert[1]/*-1*/].y;
        triangle[8] = vert[face[i+1].vert[1]/*-1*/].z;
        Vector tNorm = mat_computeNormal(triangle).normalized();
        fNorm[i+1].x = tNorm.x;  fNorm[i+1].y = tNorm.y;  fNorm[i+1].z = tNorm.z;
      }
    }
  }
  fd.close();
  // zero selected vertices
  selectedVertices.resize(nVert, 0);

  /* DEBUG
  for(int j=1;j<=nVert;j++)
	  printf("v%d: %g %g %g\n",j,vert[j].x,vert[j].y,vert[j].z);
  for(int j=1;j<=nFace;j++){
  	  printf("f%d: ",j);
  	  for(int h=0;h<face[j].nVert; h++)
  	          printf("%d ",face[j].vertId->at(h));
  	  printf("\n");
  }
  */

  /***da rimuovere ***
     for(int k=1;k<=nFace;k++){
   	  for(int i=0;i<face[k].nVert;i++){
   		  face[k].vert[i]=face[k].vertId->at(i);
   	  }
     }
  ******************************/


  this->mod_LoadSubdivMesh();

  return true;
}

/* A partire dalla Mod_SubdivMesh "mesh" con vertici e facce, 
 riempie tutti i campi previsti per la struttura Mod_SubdivMesh
 necessari per le operazioni booleane */
void Mod_SubdivMesh::mod_LoadSubdivMesh() {

  typedef struct {
    int n;
    int *vedge;
    int *epos;
    int *f1pos;
    int *f2pos;
  } work_area;

  int nvert, nedge, nface, nval;
  Vertex *vert;
  Edge *edge;
  Face *face;
  work_area *buf;
  int *count_fe;
  int coppie[2*MAXVAL],ic,il;
  int i,j,k,l,ii,jj,kk,ij,i_st,i_en,ijc;

  nvert=this->nVert;
  nface=this->nFace;
  vert=this->vert;
  face=this->face;

  #if DEBUG_LOAD
    printf("NVertices = %d \n",nvert);
    printf(" Coord Vertices \n");
    for (k=1; k<=nvert; k++)
      printf(" %d> %f %f %f\n",k,vert[k].x,vert[k].y,vert[k].z);
    printf("NFaces = %d \n",nface);
    printf(" Indici Face \n");
    for (k=1; k<=nface; k++) {
      printf(" %d> ",k);
      for (j=0; j<face[k].nVert; j++)
        printf(" %d ",face[k].vert[j]);
      printf("tan:%d\n",face[k].tan);
     }
  #endif

  /* alloco una struttura di lavoro che vado a riempire e sfruttando la
     quale posso poi riempire agevolmente le strutture EV, EF, VE */
    buf = (work_area *)malloc((this->nVert+1)*sizeof(work_area) );
    if (!buf) ErrorMessage();
    for (k=1; k<=nvert; k++) {
      buf[k].n = 0;
      buf[k].vedge = (int *)calloc( MAXVAL,sizeof(int) );
      buf[k].epos = (int *)calloc( MAXVAL,sizeof(int) );
      buf[k].f1pos = (int *)calloc( MAXVAL,sizeof(int) );
      buf[k].f2pos = (int *)calloc( MAXVAL,sizeof(int) );
      if (buf[k].vedge==NULL || buf[k].epos==NULL || buf[k].f1pos==NULL || buf[k].f2pos==NULL)
        printf("%d Non sono riuscito ad allocare memoria\n",k);
    }

    nedge=1;
    for (k=1; k<=nface; k++) {
      for (ii=0; ii<face[k].nVert; ii++) {
        if (ii<face[k].nVert-1) {
          i_st=ii;
          i_en=ii+1;
        } else {
          i_st=face[k].nVert-1;
          i_en=0;
        }
  /* controllo che non esista già memorizzato l'Edge di vertici invertiti */
        i=face[k].vert[i_en];
        kk=buf[i].n;
        j=0;
        while (j<kk && buf[i].vedge[j]!=face[k].vert[i_st])
          j++;
  /* controllo che non esista già memorizzato l'Edge di vertici dritti */
        if (j==kk) {
          i=face[k].vert[i_st];
          kk=buf[i].n;
          j=0;
          while (j<kk && buf[i].vedge[j]!=face[k].vert[i_en])
            j++;
          if (j==kk) {
            buf[i].vedge[kk]=face[k].vert[i_en];
            buf[i].epos[kk]=nedge;
            buf[i].f1pos[kk]=k;
            nedge++;
            buf[i].n++;
          }
        } else buf[i].f2pos[j]=k;
      }
    }
    nedge--;
    this->nEdge = nedge;
  /* calcolo il numero di edge della mesh come la sommatoria di buf[k].n
     ricalcolo di nedge a posteriori come controllo
       mesh->nedge = nedge = 0;
       for (k=1; k<=nvert; k++)
         nedge += buf[k].n;
     alloco la struttura Edge */
    this->edge = edge = (Edge*)malloc((nedge+1)*sizeof(Edge));
    if (!this->edge) ErrorMessage();
    if(edge == NULL)
      printf("NON sono riuscito ad allocare memoria\n");
  #if DEBUG_LOAD
    printf("NEdges = %d \n",nedge);
  #endif

  /* alloco un array di contatori temporaneo */
    count_fe = (int *)calloc((this->nVert+1), sizeof(int) );
    for (k=1; k<=nface; k++) {
      for (ii=0; ii<face[k].nVert; ii++) {
        kk=face[k].vert[ii];
        if(buf[kk].vedge[count_fe[kk]]==face[k].vert[(ii+1)%face[k].nVert]) {
          face[k].edge[ii]=buf[kk].epos[count_fe[kk]];
          count_fe[kk]++;
        } else {
          i=0;
          kk=face[k].vert[(ii+1)%face[k].nVert];
          while(i<buf[kk].n && buf[kk].vedge[i]!=face[k].vert[ii])
            i++;
            face[k].edge[ii]=buf[kk].epos[i];
        }

      }
    }
  /* Inizializzo le valenze dei vertici a zero */
    for (k=1; k<=nvert; k++) {
      vert[k].val=0;
      vert[k].bound=false;
    }
  /* Scorro la struttura di lavoro buf e riempio le strutture
    edge[].vert[] (EV)
    edge[].neigh[] (EF)
    vert[].edge[] (VE)
    face[].edge[] (FE)
    e i campi
    edge[].val (che mi informa se un edge è di bordo o meno)
    vert[].val (valenza di un vertice)
    Nota: in VE ci sono indici positivi e negativi di edge;
    se vert[i].edge[k] è positivo vuol dire che l'edge memorizzato
    ha come primo estremo il vertice di indice 'i', se invece è negativo vuol dire che
    il vertice di indice 'i' e' il secondo estremo dell'edge.
    praticamente riesco a memorizzare l'orientamento dell'edge
    che mi sembra molto importante e viene sfruttato poi per mettere
    in senso antiorario sia la lista VV che FF. */
    for (k=1; k<=nvert; k++)
      for (i=0; i<buf[k].n; i++) {
        edge[buf[k].epos[i]].vert[0]=k;
        edge[buf[k].epos[i]].vert[1]=buf[k].vedge[i];
        j=buf[k].f1pos[i];
        edge[buf[k].epos[i]].neigh[0]=j;
        edge[buf[k].epos[i]].val=1;
  /*      face[j].edge[count_fe[j]]=buf[k].epos[i]; */
  /*      count_fe[j]++; */
        j=buf[k].f2pos[i];
      if (j>0) {
        edge[buf[k].epos[i]].neigh[1]=j;
        edge[buf[k].epos[i]].val++;
  /*      face[j].edge[count_fe[j]]=buf[k].epos[i]; */
  /*      count_fe[j]++; */
      } else edge[buf[k].epos[i]].neigh[1]=0;
      vert[k].edge[vert[k].val]=buf[k].epos[i];
      vert[k].val++;
      j=buf[k].vedge[i];
      vert[j].edge[vert[j].val]=-buf[k].epos[i];
      vert[j].val++;
    }
  /* dealloco l'array count_fe */
     free(count_fe);

  /* stampa delle liste EV, EF,VE e FE */
  #if DEBUG_LOAD
    printf("\n");
    for (k=1; k<=nedge; k++) {
      printf("Edge %d: EdgeVert %d %d EdgeFace %d %d \n",k,edge[k].vert[0],edge[k].vert[1],
      edge[k].neigh[0],edge[k].neigh[1]);
    }
    printf("\n");
    for(i=1; i<=nface; i++) {
      printf("Face %d FaceEdge: ",i);
      for(j=1; j<=face[i].nVert; j++)
        printf(" %d",face[i].edge[j-1]);
        printf("\n");
    }
    printf("\n");
    for(i=1; i<=nvert; i++) {
      printf("Vertex %d VertEdges:",i);
      for (k=0; k<vert[i].val; k++)
        printf("%d ",vert[i].edge[k]);
      printf("\n");
    }
  #endif
  /* Controllo se un vertice è di bordo:
     un vertice è di bordo se appartiene a un edge di bordo, cioè con val=1
     Riempio vert.bound
     Riempio vert.val
     GC: secondo me basta fare un ciclo sugli edge;
     se l'edege i è di bordo, mettere vert[i].bound a TRUE */
    for(k=1; k<=nedge; k++)
      if (edge[k].val==1) {
        vert[edge[k].vert[0]].bound = true;
        vert[edge[k].vert[1]].bound = true;
    }

  /* stampa per ogni vertice se è di bordo o meno */
  #if DEBUG_LOAD
    printf("\n");
    for(i=1; i<=nvert; i++)
      printf("Vertice %d, Val= %d, Boundary: %d\n",i,vert[i].val,vert[i].bound);
  #endif

  /* Ora che ho la lista, per ogni vertice, degli edge (VE) li vorrei ordinare
     in senso antiorario;
     l'unico modo che vedo è costruire una catena ordinate di facce usando la lista EF
     copio gli EF in un array di coppie di indici di facce e ordino le coppie per fare
     una catena continua e a giro; tengo conto di eventuali edge di bordo controllando
     se hanno una sola faccia adiacente */
    for(i=1; i<=nvert; i++) {
      ic=0;
      for (k=0; k<vert[i].val; k++) {
        j=vert[i].edge[k];
        if (j>0) {
          if (edge[j].neigh[1]==0) {
            coppie[ic]=coppie[0];
            ic++;
            coppie[ic]=coppie[1];
            ic++;
            coppie[0]=edge[j].neigh[1];
            coppie[1]=edge[j].neigh[0];
            swap_int(&(vert[i].edge[k]),&(vert[i].edge[0]));
          } else {
            coppie[ic]=edge[j].neigh[1];
            ic++;
            coppie[ic]=edge[j].neigh[0];
            ic++;
          }
        } else {
            if (edge[-j].neigh[0]==0) {
              coppie[ic]=coppie[0];
              ic++;
              coppie[ic]=coppie[1];
              ic++;
              coppie[0]=edge[-j].neigh[0];
              coppie[1]=edge[-j].neigh[1];
              swap_int(&(vert[i].edge[k]),&(vert[i].edge[0]));
            } else {
              coppie[ic]=edge[-j].neigh[0];
              ic++;
              coppie[ic]=edge[-j].neigh[1];
              ic++;
            }
          }
        }
        for (k=0; k<vert[i].val-1; k++)
          for (j=k+1; j<vert[i].val; j++) {
            if (coppie[2*k+1]==coppie[2*j]) {
              swap_int(&coppie[2*(k+1)],&coppie[2*j]);
              swap_int(&coppie[2*(k+1)+1],&coppie[2*j+1]);
              swap_int(&(vert[i].edge[k+1]),&(vert[i].edge[j]));
              j=vert[i].val;
            }
          }

          if (vert[i].bound) {
            vert[i].face[0]=coppie[1];
            j=1;
          } else {
            vert[i].face[0]=coppie[0];
            vert[i].face[1]=coppie[1];
            j=2;
          }
          for (k=1; k<vert[i].val-1; k++) {
            vert[i].face[j]=coppie[2*k+1];
            j++;
          }

  /* A partire dai VF e dalle FV che sono entrambe gia' in senso antiorario
     estraggo nello stesso senso i neigh dei vertici (VV)
     nota che se il vertice è di bordo le facce in VF sono una in meno della valenza */
          ic=0;
          il=0;
          if (vert[i].bound)
          for (k=0; k<vert[i].val-1; k++) {
            jj=vert[i].face[k];
            if (k==0)
              ii=mod_find_in_face(face,jj,i);
            else
              ii=mod_find_in_face(face,jj,vert[i].neigh[ic-1]);
            for (l=1; l<face[jj].nVert-il; l++) {
              vert[i].neigh[ic]=face[jj].vert[(l+ii)%face[jj].nVert];
              ic++;
            }
            il=1;
          } else
            for (k=0; k<vert[i].val; k++) {
              jj=vert[i].face[k];
              if (k==0)
                ii=mod_find_in_face(face,jj,i);
              else
                ii=mod_find_in_face(face,jj,vert[i].neigh[ic-1]);
              for (l=1; l<face[jj].nVert-il; l++) {
                vert[i].neigh[ic]=face[jj].vert[(l+ii)%face[jj].nVert];
                ic++;
              }
              il=1;
              if (k==vert[i].val-2)
                il=2;
            }
            vert[i].n_neigh=ic;
          }
  /* printf("qui non arriva\n"); */

  /* stampa la lista VE ordinata */
  #if DEBUG_LOAD
     printf("\n");
     for(i=1; i<=nvert; i++) {
       printf("Vertex %d VertEdge ",i);
       for (k=0; k<vert[i].val; k++)
         printf("%d ",vert[i].edge[k]);
         printf("\n");
     }
  /* stampa la lista VF */
     printf("\n");
     for(i=1; i<=nvert; i++) {
       printf("Vertex %d VertFace ",i);
       if (vert[i].bound)
         for (k=0; k<vert[i].val-1; k++)
           printf("%d ",vert[i].face[k]);
       else
         for (k=0; k<vert[i].val; k++)
           printf("%d ",vert[i].face[k]);
           printf("\n");
     }
  /* stampa la lista VV in senso antiorario */
     printf("\n");
     for(i=1; i<=nvert; i++) {
       printf("Vertex %d VertVert ",i);
       for (k=0; k<vert[i].n_neigh; k++)
         printf("%d ",vert[i].neigh[k]);
         printf("\n");
     }
  #endif

  /* Controllo se una faccia è di bordo
     Riempio face.bound
     Una faccia è di bordo se ha un vertice di bordo
     GC: anche qui dovrebbe bastare un ciclo sui vertici e se questo
     è di bordo si etichettano tutte le facce di cui lui è vertice
     come di bordo. Sappiamo già le facce adiacenti di un vertice? */
    for(i=1; i<=nface; i++) {
      face[i].bound = false;
      j=0;
      while( (j<face[i].nVert)&&(!face[i].bound) ) {
        if(vert[face[i].vert[j]].bound)
          face[i].bound = true;
        j++;
      }
    }

  /* Cerco infine per ogni faccia, le facce che la circondano in senso
     antiorario (FF).
     se la faccia è di bordo bisogna iniziare da una intorno di bordo.
     la procedura è simile a trovare i vertici intorno ad un vertice; usero'
     cioe' la lista Face per avere i vertici di una faccia e per ognuno dalla
     lista VF estraggo a giro le facce. */
    for(i=1; i<=nface; i++) {
      ic=0;
      il=0;
      ij=0;

  /* ricerca dell'ultimo vertice di bordo della face
     una volta trovato uno di bordo, il ciclo while che
     segue scorre gli eventuali vertici di bordo
     a lui adiacenti nell'ordine antiorario */

  /* Versione errata 1:
    for (j=0; j<face[i].n_v_e; j++)
        if(vert[face[i].vert[j]].bound)
          ij=j;
  */

  /* Versione errata 2:
      for (j=0; j<face[i].n_v_e; j++)
        if (vert[face[i].vert[j]].bound)
          ij=j;
      if (ij<face[i].n_v_e-1)  evita loop infiniti con mesh costituite da un unica faccia
        while(vert[face[i].vert[(ij+1)%face[i].n_v_e]].bound)
          ij=(ij+1)%face[i].n_v_e;
  */

      if (face[i].bound) {
       /* ricerca dell'ultimo vertice di bordo della face
          una volta trovato uno di bordo, il ciclo while che
          segue scorre gli eventuali vertici di bordo
          a lui adiacenti nell'ordine antiorario. */
        ij=0;
        while(!vert[face[i].vert[ij]].bound)
          ij++;

  /* GC08 aggiunto contatore ijc per evitare ciclo infinito se tutti i
     vertici di una faccia sono di bordo */
        ijc=1;
        while(vert[face[i].vert[(ij+1)%face[i].nVert]].bound &&
              ijc<=face[i].nVert) {
            ij=(ij+1)%face[i].nVert;
        ijc++;
        }
      }
  /* GC08 versione nuova per gestire ring di una faccia di bordo con
     discontinuita' delle facce intorno */
      for (j=0; j<face[i].nVert; j++) {
        jj=face[i].vert[(j+ij)%face[i].nVert];
        if (vert[jj].bound)
          nval=vert[jj].val-1;
        else
          nval=vert[jj].val;
  /* GC08 cerco l'indice della faccia attuale nel ring di facce del vertice jj */
        ii=mod_find_in_vert_face(vert,jj,i);
       /* if(i==24){
           for(int s=0;s<face[i].nVert;s++)
        	   printf("\nj: %d Con Edge %d \n",j,face[i].edge[s]);
        printf("Faccia:%d Facce a cui appartiene il vertice %d :\n",i,jj);
           for (l=1; l<nval; l++)
        		printf("face: %d\n",vert[jj].face[(l+ii)%nval]);
        }*/
  /* GC08 esamino le facce del ring del vertice jj a partire dalla
     successiva alla i e la inserisco se non gia' precedentemente inserita */
        for (l=1; l<nval; l++) {
         if (ic!=0) {
            if(vert[jj].face[(l+ii)%nval]!=face[i].neigh[ic-1]) {
              face[i].neigh[ic]=vert[jj].face[(l+ii)%nval];
              ic++;
            }
         } else {
            face[i].neigh[ic]=vert[jj].face[(l+ii)%nval];
            ic++;
          }
        }
      }

  /* GC08 test finale per evitare che l'ultima faccia sia ripetuta
     rispetto alla prima */
       if(ic>1 && face[i].neigh[0]==face[i].neigh[ic-1])
            ic--;

       face[i].val=ic;
    }

  /* stampa */
  #if DEBUG_LOAD
    printf("\n");
    for (i=1; i<=nface; i++) {
      printf("Face %d: FaceFace: ",i);
      for (ic=0; ic<face[i].val; ic++)
        printf("%d ",face[i].neigh[ic]);
        printf("\n");
    }
  #endif

    // compute face normal
    fNorm = (Vertex*) malloc((nface+1)*sizeof(Vertex));
    if (!fNorm)
    {
    	cout << "Error: Malloc fails -> exit!!" << endl;
    	exit(2);
    }
    for (i=0; i<nface; i++) {
    	float triangle[9];
    	triangle[0] = vert[face[i+1].vert[0]/*-1*/].x;
    	triangle[1] = vert[face[i+1].vert[0]/*-1*/].y;
    	triangle[2] = vert[face[i+1].vert[0]/*-1*/].z;
    	triangle[3] = vert[face[i+1].vert[2]/*-1*/].x;
    	triangle[4] = vert[face[i+1].vert[2]/*-1*/].y;
    	triangle[5] = vert[face[i+1].vert[2]/*-1*/].z;
    	triangle[6] = vert[face[i+1].vert[1]/*-1*/].x;
    	triangle[7] = vert[face[i+1].vert[1]/*-1*/].y;
    	triangle[8] = vert[face[i+1].vert[1]/*-1*/].z;
    	Vector tNorm = mat_computeNormal(triangle).normalized();
    	fNorm[i+1].x = tNorm.x;  fNorm[i+1].y = tNorm.y;  fNorm[i+1].z = tNorm.z;
    }


  /* dealloco struttura di lavoro */
    for (k=1; k<=nvert; k++) {
      free(buf[k].vedge);
      free(buf[k].epos);
      free(buf[k].f1pos);
      free(buf[k].f2pos);
    }
    free(buf);

    /* marco gli edge tangenti */
    for(i=1;i<=nface;i++){
    	for(j=0;j<face[i].nVert;j++)
    		edge[face[i].edge[j]].tan = 0;
    }
    for(i=1;i<=nface;i++){
    	if(face[i].tan==1){
    		for(j=0;j<face[i].nVert;j++){
    			edge[face[i].edge[j]].tan = 1;
    		}
    	}

    }

    /* da rimuovere*
    for(i=1;i<=nface;i++){
    	face[i].vertId = new vector<int>(face[i].nVert);
    	 for (j=0; j<face[i].nVert; j++)
    	    	face[i].vertId->at(j) = face[i].vert[j];
    }
    */
    return;
  //    return mesh;
}

// save the obj model
bool Mod_SubdivMesh::mod_save()
{
  return mod_save(filename, filepath);
}

bool Mod_SubdivMesh::mod_save(string fName, string fPath)
{
  // obj wavefront to load
  string file = "";
  file += fPath;
  file += fName;
  file += "_res.obj";

  ofstream fd;
  fd.open(file.c_str());
  // open file to write
  if (!fd.good())
  {
    cout << "Error: Unable to open " << filename << " -> exit!!" << endl;
    exit(1);
  }
  // notehead
  fd << "#\t  Wavefront obj file" << endl;
  fd << "#\tCreated using Craftmesh" << endl;
  fd << endl;
  fd << "#\t University of Bologna" << endl;

  // vertices
  for (int i=1; i<=nVert; i++)
    fd << "v " << vert[i].x << " " << vert[i].y << " " << vert[i].z << endl;

  // faces
  for (int i=1; i <= nFace; i++)
  {
    fd << "f";
    for (int j=0; j<face[i].nVert; j++)
      fd << " " << face[i].vert[j];//Id->at(j);
    fd << endl;
  }

  fd.close();

  return true;
}

// copy source mesh in another mesh object
void Mod_SubdivMesh::mod_copyMesh(Mod_SubdivMesh *Mdst) {

  int i, j;

#if DEBUG_GLUE
  printf("COPY---------------------------->>\n");
#endif

  /* copia della mesh sorgente in mesh destinazione */
  Mdst->nVert = this->nVert;
  Mdst->nFace = this->nFace;
  Mdst->vert  = (Vertex*) malloc((Mdst->nVert + 1) * sizeof(Vertex));
  if (!Mdst->vert) ErrorMessage();
  Mdst->face = (Face*) malloc((Mdst->nFace + 1) * sizeof(Face));
  if (!Mdst->face) ErrorMessage();

  /* ... copia dei vertici e ... */
  for (i=1; i <= Mdst->nVert; i++) {
    Mdst->vert[i].x = this->vert[i].x;
    Mdst->vert[i].y = this->vert[i].y;
    Mdst->vert[i].z = this->vert[i].z;
    Mdst->vert[i].tan = this->vert[i].tan;
    Mdst->vert[i].ins = this->vert[i].ins;
#if DEBUG_GLUE
    if (i==1) printf("Lista vertici:\n");
    printf("%f %f %f\n",Mdst->vert[i].x,Mdst->vert[i].y,Mdst->vert[i].z);
#endif
  }

  /* ... copia delle faccette. */
  for (i=1; i <= Mdst->nFace; i++) {

    Mdst->face[i].nVert = this->face[i].nVert;
    Mdst->face[i].pointInto = this->face[i].pointInto;
    Mdst->face[i].tan = this->face[i].tan;
    Mdst->face[i].bb = this->face[i].bb;

#if DEBUG_GLUE
    if (i==1) printf("Lista facce:\n");
    printf("Face no.%d NVertices  %d \n",i,Mdst->face[i].nVert);
#endif
    for (j=0; j < Mdst->face[i].nVert; j++)
      Mdst->face[i].vert[j] = this->face[i].vert[j];
  }

  Mdst->mod_LoadSubdivMesh();

#if DEBUG_GLUE
  printf("<<----------------------------COPY\n");
#endif

  return;
}

// set the mesh center coordinates
void Mod_SubdivMesh::mod_setCenterAndRadius(void)
{
  // set the center of the mesh
  float bbMinX, bbMaxX,
        bbMinY, bbMaxY,
        bbMinZ, bbMaxZ;
  bbMinX = bbMinY = bbMinZ = AMAXVALUE;  // numeric_limits<float>::max();
  bbMaxX = bbMaxY = bbMaxZ = AMINVALUE;  // numeric_limits<float>::min();
  for (int i=1; i<=nVert; i++)
  {
    if (vert[i].x < bbMinX) bbMinX = vert[i].x;
    if (vert[i].x > bbMaxX) bbMaxX = vert[i].x;
    if (vert[i].y < bbMinY) bbMinY = vert[i].y;
    if (vert[i].y > bbMaxY) bbMaxY = vert[i].y;
    if (vert[i].z < bbMinZ) bbMinZ = vert[i].z;
    if (vert[i].z > bbMaxZ) bbMaxZ = vert[i].z;
  }
  center.x = (bbMaxX-bbMinX)/2.0 + bbMinX;
  center.y = (bbMaxY-bbMinY)/2.0 + bbMinY;
  center.z = (bbMaxZ-bbMinZ)/2.0 + bbMinZ;
  // set the sphere radius to capture the mesh
  float minD, maxD;
  Vector minDistance(bbMinX, bbMinY, bbMinZ);
  Vector maxDistance(bbMaxX, bbMaxY, bbMaxZ);
  minD = center.distanceTo(minDistance);
  maxD = center.distanceTo(maxDistance);
  rCapture = (maxD + minD)/2.0;
  if (rCapture < MIN_RCAPTURE) rCapture = MIN_RCAPTURE;
  else if (rCapture > MAX_RCAPTURE) rCapture = MAX_RCAPTURE;

  return;
}

// update the mesh center coordinates
void Mod_SubdivMesh::mod_translateCenter(Mat_Vector trans)
{
  mat_translateVector(center, trans, center);
  return;
}

// draw the model
void Mod_SubdivMesh::mod_drawMesh(void)
{
  glDisable(GL_LIGHTING);
  glColor3fv(GL_BLU_COLOR);
  glLineWidth(2.0);
  glTranslatef(center.x, center.y, center.z);
  glMultMatrixf(rotat);
  glMultMatrixf(scale);
  glTranslatef(-center.x, -center.y, -center.z);
  glTranslatef(transl.x, transl.y, transl.z);

  // edges
  for (int i=1; i<=nFace; i++)
  {
    for (int j=0; j<face[i].nVert; j++)
    {
      int k = (j+1)%face[i].nVert;
      glBegin(GL_LINES);
        glVertex3f(vert[face[i].vert[j]/*-1*/].x,
                   vert[face[i].vert[j]/*-1*/].y,
                   vert[face[i].vert[j]/*-1*/].z);
        glVertex3f(vert[face[i].vert[k]/*-1*/].x,
                   vert[face[i].vert[k]/*-1*/].y,
                   vert[face[i].vert[k]/*-1*/].z);
      glEnd();
    }
  }
  glDisable(GL_DEPTH_TEST);
  glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
  glEnable(GL_BLEND);
  glColor4fv(this->color);
  // faces
  for (int i=1; i<=nFace; i++)
  {
    for (int j=0; j<face[i].nVert-2; j++)  // NOTE: run-time triangulation, it fails with concave faces
    {
      glBegin(GL_TRIANGLES);
        glVertex3f(vert[face[i].vert[0]/*-1*/].x,
                   vert[face[i].vert[0]/*-1*/].y,
                   vert[face[i].vert[0]/*-1*/].z);
        glVertex3f(vert[face[i].vert[j+1]/*-1*/].x,
                   vert[face[i].vert[j+1]/*-1*/].y,
                   vert[face[i].vert[j+1]/*-1*/].z);
        glVertex3f(vert[face[i].vert[j+2]/*-1*/].x,
                   vert[face[i].vert[j+2]/*-1*/].y,
                   vert[face[i].vert[j+2]/*-1*/].z);
      glEnd();
    }
  }
  glTranslatef(-transl.x, -transl.y, -transl.z);
  glTranslatef(center.x, center.y, center.z);
  glMultMatrixf(invsca);
  glMultMatrixf(invrot);
  glTranslatef(-center.x, -center.y, -center.z);
  glDisable(GL_BLEND);
  glEnable(GL_DEPTH_TEST);
  glEnable(GL_LIGHTING);
  return;
}

// draw the model for anaglyph vision
void Mod_SubdivMesh::mod_drawMeshSolid(void)
{
  glColor3fv(GL_GREY_COLOR);
  glTranslatef(center.x, center.y, center.z);
  glMultMatrixf(rotat);
  glMultMatrixf(scale);
  glTranslatef(-center.x, -center.y, -center.z);
  glTranslatef(transl.x, transl.y, transl.z);

  // faces
  for (int i=1; i<=nFace; i++)
  {
    for (int j=0; j<face[i].nVert-2; j++)  // NOTE: run-time triangulation, it fails with concave faces
    {
      glNormal3f(fNorm[i].x, fNorm[i].y, fNorm[i].z);
      glBegin(GL_TRIANGLES);
        glVertex3f(vert[face[i].vert[0]/*-1*/].x,
                   vert[face[i].vert[0]/*-1*/].y,
                   vert[face[i].vert[0]/*-1*/].z);
        glVertex3f(vert[face[i].vert[j+1]/*-1*/].x,
                   vert[face[i].vert[j+1]/*-1*/].y,
                   vert[face[i].vert[j+1]/*-1*/].z);
        glVertex3f(vert[face[i].vert[j+2]/*-1*/].x,
                   vert[face[i].vert[j+2]/*-1*/].y,
                   vert[face[i].vert[j+2]/*-1*/].z);
      glEnd();
    }
  }

  glTranslatef(-transl.x, -transl.y, -transl.z);
  glTranslatef(center.x, center.y, center.z);
  glMultMatrixf(invsca);
  glMultMatrixf(invrot);
  glTranslatef(-center.x, -center.y, -center.z);
  return;
}

// draw the catching sphere of the mesh
void Mod_SubdivMesh::mod_drawCatchingSphere(void)
{
  glDisable(GL_LIGHTING);
  glColor3fv(GL_GREY_COLOR);
  glLineWidth(0.1);
  glTranslatef(center.x, center.y, center.z);
  glRotatef(90.0, 1.0, 0.0, 0.0);
  glutWireSphere(rCapture, 16, 8);
  glRotatef(-90.0, 1.0, 0.0, 0.0);
  glTranslatef(-center.x, -center.y, -center.z);
  glEnable(GL_LIGHTING);
  return;
}

// draw the reference cylinder of the mesh
void Mod_SubdivMesh::mod_drawRefCylinder(int visionMode)
{
  glDisable(GL_LIGHTING);
  glDisable(GL_DEPTH_TEST);
  glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
  glEnable(GL_BLEND);
  if (visionMode == GL_VMONO || visionMode == GL_VMULTI) glColor4fv(this->color);
  else glColor4fv(GL_GREY_CYLI);
  // on XZ plane
  glTranslatef(center.x, 0.0, center.z);
  glRotatef(90.0, 1.0f, 0.0f, 0.0f);
  glutSolidCylinder(rCapture, 0.2, 32, 32);
  glRotatef(-90.0, 1.0f, 0.0f, 0.0f);
  glTranslatef(-center.x, 0.0, -center.z);
  // on XY plane
  glTranslatef(center.x, center.y, -150.0);
  glutSolidCylinder(rCapture, 0.2, 32, 32);
  glTranslatef(-center.x, -center.y, 150.0);
  glDisable(GL_BLEND);
  glEnable(GL_DEPTH_TEST);
  glEnable(GL_LIGHTING);
  return;
}

// apply the transformations
void Mod_SubdivMesh::mod_applyTransformations(void)
{
  // apply transformations (S)cale (R)otation (T)ranslation [(v+T-c')*S*R+c']
  float vertexTemp[3];
  for (int i=1; i<=nVert; i++)
  {
    vertexTemp[0] = vert[i].x + transl.x - center.x;
    vertexTemp[1] = vert[i].y + transl.y - center.y;
    vertexTemp[2] = vert[i].z + transl.z - center.z;
    mat_vectorMatrixMultiply(vertexTemp, scale, vertexTemp);
    mat_vectorMatrixMultiply(vertexTemp, rotat, vertexTemp);
    vert[i].x = vertexTemp[0] + center.x;
    vert[i].y = vertexTemp[1] + center.y;
    vert[i].z = vertexTemp[2] + center.z;
  }

  // reset transformations
  transl.x = 0.0;  transl.y = 0.0;  transl.z = 0.0;
  mat_glMatrixfIdentity(scale);
  mat_glMatrixfIdentity(invsca);
  mat_glMatrixfIdentity(rotat);
  mat_glMatrixfIdentity(invrot);

  return;
}

// draw the selected vertices
void Mod_SubdivMesh::mod_drawSelectedVertices(int visionMode)
{
  glDisable(GL_LIGHTING);
  glTranslatef(center.x, center.y, center.z);
  glMultMatrixf(rotat);
  glMultMatrixf(scale);
  glTranslatef(-center.x, -center.y, -center.z);
  glTranslatef(transl.x, transl.y, transl.z);

  if (visionMode == GL_VMONO || visionMode == GL_VMULTI) glColor3fv(GL_YELLOW_COLOR);
  else glColor3fv(GL_GREY_COLOR);

  for (unsigned int i=0; i<selectedVertices.size(); i++)
  {
    if (selectedVertices[i])
    {
      //i-->i+1
      glTranslatef(vert[i+1].x, vert[i+1].y, vert[i+1].z);
      glutSolidSphere(1, 5, 5);
      glTranslatef(-vert[i+1].x, -vert[i+1].y, -vert[i+1].z);
    }
  }

  glTranslatef(-transl.x, -transl.y, -transl.z);
  glTranslatef(center.x, center.y, center.z);
  glMultMatrixf(invsca);
  glMultMatrixf(invrot);
  glTranslatef(-center.x, -center.y, -center.z);
  glEnable(GL_LIGHTING);
  return;
}

// search in a face for a given vertex index
int Mod_SubdivMesh::mod_find_in_face(Face *face,int ind_face, int ind_vert) {

  int k;
  for (k=0; k<face[ind_face].nVert; k++)
    if (face[ind_face].vert[k]==ind_vert)
    break;
  return(k);
}

// search a vertex with a given face index
int Mod_SubdivMesh::mod_find_in_vert_face(Vertex *vert,int ind_vert, int ind_face) {

  int k,nval;
  if (vert[ind_vert].bound)
    nval=vert[ind_vert].val-1;
  else
    nval=vert[ind_vert].val;
  for (k=0; k<nval; k++)
    if (vert[ind_vert].face[k]==ind_face)
    break;
  return(k);
}

// search a vertex index with a given vertex
int Mod_SubdivMesh::mod_find_in_vert(Vertex v){

	int res = -1;

	for(int i=1;i<=nVert;i++){
		if(vert[i].x==v.x && vert[i].y==v.y && vert[i].z==v.z){
			res = i;
			break;
		}
	}

	return res;
}

// Funzioni per l'interrogazione degli elementi costitutivi di una mesh
void Mod_SubdivMesh::mod_inquiryV(int index) {

  int i;

  printf("========================================\n");
  printf("Le coordinate del vertice indicato sono:\n");
  printf("X: %f\tY: %f\tZ: %f\n",
         vert[index].x,vert[index].y,vert[index].z);

  if (vert[index].bound) printf("è un vertice di bordo e ha valenza");
  else printf("non è un vertice di bordo e ha valenza");
  printf(" %d;\n", vert[index].val);

  printf("su di esso incidono i seguenti edge:");
  for (i=0; i < MAXVAL; i++) {
    if (!vert[index].edge[i]) break;
    else printf(" %d",vert[index].edge[i]);
  }
  printf("\nil - indica edge entrante;\n");

  printf("tale vertice è costitutivo delle seguenti face:");
  for (i=0; i < MAXVAL; i++) {
    if (!vert[index].face[i]) break;
    else printf(" %d",vert[index].face[i]);
  }
  printf(";\n");

  printf("i vertici del primo ring attorno ad esso sono %d:",
          vert[index].n_neigh);
  for (i=0; i < MAXVAL; i++) {
    if (!vert[index].neigh[i]) break;
    else printf(" %d",vert[index].neigh[i]);
  }
  printf(".\n========================================\n");

  return;
}

void Mod_SubdivMesh::mod_inquiryE(int index) {

  printf("========================================\n");
  printf("L'edge indicato è definito da un vertex\n");
  printf("uscente e da uno entrante: %d %d\n",
          edge[index].vert[0],edge[index].vert[1]);

  if (!edge[index].neigh[1]) printf("è un edge di bordo e ha valenza");
  else printf("non è un edge di bordo e ha valenza");
  printf(" %d;\n", edge[index].val);

  printf("tale edge è costitutivo delle seguenti face: %d",
          edge[index].neigh[0]);
  if (edge[index].neigh[1]) printf(" %d.\n",edge[index].neigh[1]);
  else printf(".\n");

  printf("========================================\n");
  return;
}

void Mod_SubdivMesh::mod_inquiryF(int index) {

  int i;
  int *v;
  bool t;
  double param[4];  /* coefficienti del piano per il test di complanarità */

  printf("========================================\n");
  printf("La face indicata è definita da %d vertex:",face[index].nVert);
  for (i=0; i < MAXVAL; i++) {
    if (!face[index].vert[i]) break;
    else printf(" %d",face[index].vert[i]);
  }
  if (gen_complaneTest(this, index, param))
    printf("\ntutti complanari fra loro;\n");
  else printf("\nnon complanari fra loro;\n");

  printf("conseguentemente dai seguenti %d edge:",face[index].nVert);
  for (i=0; i < MAXVAL; i++) {
    if (!face[index].edge[i]) break;
    else printf(" %d",face[index].edge[i]);
  }
  printf(";\n");

  /* ATTENZIONE: il test che individua eventuali incroci tra edge potrebbe
                 fallire se la face è concava. */
  if (gen_crossEdge(this, index)) printf("rispetta le regole Two-Manifold, ed è ");
  else printf("presenta un incrocio tra gli edge, ed è ");

  v = gen_shapeFace(this, index, &t);
  if (t) printf("convessa;\n");
  else {
    printf("concava a causa del vertice");
    for (i=0; i < face[index].nVert; i++)
      if (v[i])
        printf(" %d",v[i]);
    printf(";\n");
  }
  free(v);

  if (face[index].bound) printf("è una face di bordo e ha valenza");
  else printf("non è un face di bordo e ha valenza");
  printf(" %d;\n", face[index].val);

  printf("le face del primo ring sono:");
  for (i=0; i < MAXVAL; i++) {
    if (!face[index].neigh[i]) break;
    else printf(" %d",face[index].neigh[i]);
  }

  printf(".\n========================================\n");

  return;
}

// Print an error message
void Mod_SubdivMesh::ErrorMessage() {
      printf("ERROR: NO memory!!\n");
      exit(1);
}

// swap two int
void Mod_SubdivMesh::swap_int(int *a, int *b) {

  int temp;

  temp=*a;
  *a=*b;
  *b=temp;
}
